Maine Coastal Current RSTARS

Detailing the Regime Shifts in Maine Coastal Current Behavior

Published

March 14, 2025

Code 
{
  library(sf) 
  library(fvcom) 
  library(tidyverse)
  library(gmRi)
  library(patchwork)
  library(rnaturalearth)
  library(showtext)
  library(ncdf4)
  # Cyclic color palettes in scico
  # From: https://www.fabiocrameri.ch/colourmaps/
  library(scico)
  library(legendry)
}

# namespace conflicts
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")


# Set the theme
theme_set(theme_classic() + map_theme())

# Project paths
lob_ecol_path <- cs_path("mills", "Projects/Lobster ECOL")
fvcom_path    <- cs_path("res", "FVCOM/Lobster-ECOL")
poly_paths    <- cs_path("mills", "Projects/Lobster ECOL/Spatial_Defs")

# Shapefiles
new_england <- ne_states("united states of america", returnclass = "sf") %>% 
  filter(postal %in% c("VT", "ME", "RI", "MA", "CT", "NH", "NY", "MD", "VA", "NJ", "DE", "NC", "PA", "WV"))
canada <- ne_states("canada", returnclass = "sf")

# # # Support functions for FVCOM
source(here::here("R/FVCOM_Support.R"))
Code 
# Use GMRI style
use_gmri_style_rmd()

Maine Coastal Current Regime Shift Evaluation

The Maine Coastal Current can be split into a western branch (the western maine coastal current WMCC) & an eastern branch (the eastern Maine coastal current). The two branches are separated by an area off the coast of Penobscot Bay. This is a junction where a portion of the MCC will at times deflect to the South, a seasonal behavior that is influenced by river discharge.

Code 
# Project path
fvcom_processed_path <- cs_path("mills", "Projects/Lobster ECOL/FVCOM_processed/MaineCoastalCurrent")

# Load the daily surface currents (u,v)
daily_mcc_sc <- read_csv(here::here("local_data/MCC_PC_timeseries/MCC_daily_surface_velocities.csv"))

# Load the PCA Timeseries
mcc_pc_monthly <- read_csv(here::here("local_data/MCC_PC_timeseries/monthly_surface_velocities_PCA_timeseries.csv"))
Code 
# Loading FVCOM


# Get some file names
fvcom_surfbot_files <- setNames(
  list.files(fvcom_path, full.names = T, pattern = ".nc"),
  str_remove(list.files(fvcom_path, full.names = F, pattern = ".nc"), ".nc"))

# Open one
gom3_early <- nc_open(fvcom_surfbot_files[[1]])

# Get the mesh itself as a simple feature collection
gom3_mesh <- get_mesh_geometry(gom3_early, what = 'lonlat')
Code 
#Original Area
mcc_turnoff_coords <- tribble(
  ~"lon", ~"lat",
  -69.34, 43.21,   # Bottom Left
  -69.78, 43.62,   # Off Popham
  -67.45, 44.34,   # Off Jonesport
  -67.4, 43.8,     # Bottom right
  -69.34, 43.21,   # Bottom left again
)


# Make it a polygon
mcc_turnoff_poly <- st_polygon(
    list(cbind(mcc_turnoff_coords$lon, mcc_turnoff_coords$lat))) %>% 
  st_sfc(crs = 4326) %>% 
  st_as_sf() %>% 
  mutate(area = "Maine Coastal Current Region") 
st_geometry(mcc_turnoff_poly) <- "geometry"

# Project polygon
poly_projected <- st_transform(mcc_turnoff_poly, st_crs(gom3_mesh)) 

# Turn off s2
sf::sf_use_s2(FALSE)

# Trim it
mcc_studyarea_mesh <- mesh_trim(
  mesh = gom3_mesh, 
  domain = st_as_sf(poly_projected) )

Long-term Monthly Current Directions

Before proceeding too far, here is what the average flow direction is across the domain of interest to us:

Code 
# Quick verification that directions make sense
average_directions <-  daily_mcc_sc %>% 
  mutate(period = lubridate::month(time, label = TRUE)) %>% 
  bind_rows(mutate(daily_mcc_sc, period = "Year Round")) %>% 
  mutate(period = factor(period, levels = c(month.abb, "Year Round"))) %>% 
  group_by(period, lonc, latc) %>% 
  summarise(
    across(.cols = c(u,v), .fns = ~mean(.x, na.rm = T))) %>% 
  mutate(
    # Radian to degree conversion
    angle = atan2(v, u) * 180 / pi,   # Convert from radians to degrees
    angle = if_else(angle<0, angle+360, angle),
    speed = sqrt(u^2 + v^2),          # Compute speed (magnitude)
    dx = u / speed * 0.05,            # Scale x-component for visualization
    dy = v / speed * 0.05)            # Scale y-component for visualization



# Map them
yr_round_map <- average_directions %>% 
  filter(period == "Year Round") %>% 
  ggplot() +
  geom_sf(data = new_england) +
  geom_segment(
    aes(lonc, latc, xend = lonc + dx, yend = latc + dy, color = angle),
    arrow = arrow(length = unit(0.075, "cm")),
    linewidth = 0.5) +
  coord_sf(
    xlim = c(-70, -67.3), 
    ylim = c(43.2, 44.7), 
    expand = F) +
  scale_y_continuous(breaks = seq(43, 45, .5)) +
  scale_x_continuous(breaks = seq(-71, -67, 1)) +
  facet_wrap(~period, ncol = 4) +
  scale_color_scico(
    palette = 'romaO',
    breaks = c(45, 90, 135, 180, 225, 270, 315, 360),
    labels = c("NE", "N", "NW", "W", "SW", "S", "SE", "E"),
    limits = c(0,360)) +
  guides(color = guide_colring(
    reverse = T, 
    start = pi/2)) +
  theme(
    legend.position = "bottom", 
    legend.title.position = "top") +
  labs(
    x = "Longitude", y = "Latitude",
    color = "Current Flow")


yr_round_map

Code 
# Map the months independently
monthly_flow <- average_directions %>% 
  filter(period != "Year Round") %>% 
  ggplot() +
  geom_sf(data = new_england) +
  geom_segment(
    aes(lonc, latc, xend = lonc + dx, yend = latc + dy, color = angle),
    arrow = arrow(length = unit(0.03, "cm")),
    linewidth = 0.25) +
  coord_sf(
    xlim = c(-70, -67.3), 
    ylim = c(43.2, 44.7), 
    expand = F) +
  scale_y_continuous(breaks = seq(43, 45, .5)) +
  scale_x_continuous(breaks = seq(-71, -67, 1)) +
  facet_wrap(~period, ncol = 3) +
  scale_color_scico(
    palette = 'romaO',
    breaks = c(45, 90, 135, 180, 225, 270, 315, 360),
    labels = c("NE", "N", "NW", "W", "SW", "S", "SE", "E"),
    limits = c(0,360)) +
  guides(color = guide_colring(
    reverse = T, 
    start = pi/2)) +
  theme(legend.position = "none") +
  labs(
    x = "Longitude", y = "Latitude",
    color = "Current Flow")


monthly_flow

Review the MCC Principal Components

The two principal components we’re using are derived from daily FVCOM surface-layer current vector information. These first two components explain the following proportions of the variance in the current directions:

Code 
# Load the PCA we ran:
mcc_pca <- readRDS(here::here("local_data/MCC_PC_timeseries/daily_surface_velocities_PCA.rds"))

# Summary - Proportion of variance
mcc_pca_summ <- summary(mcc_pca)
mcc_pca_summ$importance[1:2, 1:2] %>% 
  round(2) %>% 
  as.data.frame() %>% 
  rownames_to_column("Feature")  %>% 
  gt::gt()
Feature PC1 PC2
Standard deviation 21.46 18.97
Proportion of Variance 0.39 0.31

The loadings for each principal component show these patterns when plotted in space:

Code 
# Pull + Reshape the Rotations / Loadings
mcc_loadings <- data.frame(
  "PC1" = mcc_pca$rotation[,"PC1"],
  "PC2" = mcc_pca$rotation[,"PC2"]) %>% 
  rownames_to_column("loc") %>% 
  separate(col = "loc", into = c("var", "elem"), sep = "_") %>% 
  pivot_longer(
    cols = starts_with("PC"), 
    names_to = "PC", 
    values_to = "PC_rotation") %>% 
  mutate(longname = if_else(var == "u", "Eastward Velocity (u)", "Northward Velocity (v)"))



# Map The Loadings
loadings_map <- mcc_studyarea_mesh %>% 
  mutate(elem = as.character(elem)) %>% 
  left_join(mcc_loadings) %>% 
  st_as_sf() %>% 
  ggplot() +
  geom_sf(data = new_england) +
  geom_sf(aes(fill = PC_rotation)) +
  facet_grid(PC~longname, labeller = label_wrap_gen(width = 10)) +
  coord_sf(
    xlim = c(-70.2, -67.1), 
    ylim = c(43.2, 44.7), 
    expand = F) +
  scale_y_continuous(breaks = seq(43, 45, .5)) +
  scale_x_continuous(breaks = seq(-71, -67, 1)) +
  scale_fill_distiller(palette = "RdBu", limits = c(-0.06, 0.06)) +
  map_theme(
    legend.title.position = "top",
    legend.title = element_text(hjust = 0.5),
    legend.position = "right",
    legend.direction = "horizontal") +
  labs(fill = "Principal Component Rotation",
       title = "MCC Principal Component Loadings/Rotations")





# Plot daily with the monthly smooth
mcc_pc_ts <- mcc_pc_monthly %>% 
  ggplot() +
  geom_line(
    aes(date, vals, color = `Principal Component`),
    linewidth = 0.4, alpha = 0.4) +
  geom_line(
    aes(date, vals_rollmean, color = `Principal Component`),
    linetype = 1) +
  rcartocolor::scale_color_carto_d() + 
  labs(
    title = "MCC Surface Current PCA Timeseries",
    y = "Principal Component Value",
    x = "Date",
    color = "12-month average")



loadings_map / mcc_pc_ts

We can see from these maps that each mesh element contributed a similar amount to each principal component, which East/West velocities contributing more to PC1 & North/South contributing more to PC2.

Further investigation into how these relate to the southward turnoff flow can be found in MCC_Daily_Workup.Qmd

STARS Regime Shifts - {rstars} Implementation

To investigate potential regime shift dynamics, each monthly principal component timeseries will be evaluated using the STARS methodology.

The following results extend the detailed approach of Stirnimann et al. 2019’s {rstars} repository.

Code 
# Stirnimann used these values in their paper:
# l = 5, 10, 15, 17.5 years, with monthly data
# Huber = 1
# Subsampling = (l + 1) / 3


# Load the function(s)
source(here::here("rstars-master","rSTARS.R"))

Regime Shifts in Maine Coastal Current Principal Components

Code 
# Get the for the
mcc_rstars <- mcc_pc_monthly %>% 
  select(`Principal Component`, time = date, vals) %>% 
  split(.$`Principal Component`) %>% 
  map_dfr(function(x){
    
    rstars(
      data.timeseries = as.data.frame(
        x[,c("time", "vals")]),
      l.cutoff = 12*7, # Seven years
      pValue = 0.05,
      Huber = 1,
      Endfunction = T,
      preWhitening = T,
      OLS = F,
      MPK = T,
      IP4 = F,
      SubsampleSize = (12*7 + 1)/3,
      returnResults = T)  },
    .id = "PC"
  ) %>% 
  mutate(
    var = "Maine Coastal Current",
    PC = str_replace_all(PC, "PC", "Principal Component "),
    shift_direction = if_else(RSI>0, "Shift Up", "Shift Down"))  %>% 
  group_by(PC) %>% 
  arrange(time) %>% 
  mutate(vals_rollmean = zoo::rollmean(vals, k = 12, fill = NA, align = "center")) %>% 
  ungroup()

Summary Plot

Code 
# Summarise the breakpoint locations
shift_points <- mcc_rstars %>% 
  filter(RSI != 0) %>% 
  dplyr::select(time, PC, var, shift_direction)


# Plot the breaks over the monthly data
ggplot() +
   geom_vline(
    data = shift_points,
    aes(
      xintercept = time,
      color = shift_direction),
    linewidth = 1.5) +
  geom_line(
    data = mcc_rstars,
     aes(time, vals, group = PC),
    linewidth = 0.2, alpha = 0.5) +
  geom_line(
    data = mcc_rstars,
     aes(time, vals_rollmean, group = PC),
    linewidth = 0.8, alpha = 0.75) +
 
  scale_color_gmri() +
  facet_wrap(~var*PC, ncol = 1, labeller = label_wrap_gen()) +
  theme(
    legend.position = "bottom",
    strip.text.y = element_text(angle = 0)) + 
  labs(color = "",
       y = "Measurement",
       title = "STARS Regime Shifts, Monthly Maine Coastal Current")

Export

Code 
# # Save it
write_csv(
  mcc_rstars,
  here::here("rstars_results/mcc_monthly_shifts.csv"))